%% Estimate stochastic part of income process (Table 2 in paper)
% Annika Bacher
% following lecture notes of G. Violante
% June 2024

%% housekeeping
clear;
clc;
close all;
tic

% !! ADD PATH TO WHERE EMPIRICAL RESULTS ARE STORED !!
datafolder = '...';


%% outer loop: seperately for men women, single and married
 gender2 = [11 12 13 14];

for gender= gender2
    if gender == 11
 csvread([datafolder '/COVmale_coup_labinc.csv']);
     elseif gender == 12
 csvread([datafolder '/COVfemale_coup_labinc.csv']);
     elseif gender == 13
 csvread([datafolder '/COVmale_sing_labinc.csv']);
     elseif gender == 14
 csvread([datafolder '/COVfemale_sing_labinc.csv']);
    end
    
data = ans;
clear ans

data(data==0) = NaN;

% prepare data to generate C matrix

P=19 ; % amount of waves available to study
T=29+1 ; % amout of calendar years the study covers + 1 for identifier

%%%%%%%%
% interpolate (linear) income in missing years (bi-annuality of PSID after 1997)
data2= zeros(length(data),T) ;
data2(:,1)=data(:,1) ;

for i=2:9
data2(:,i)=data(:,i) ;   
end 

t=0;
for i = 10:2:T
    data2(:,i)=data(:,i-t) ;
    t=t+1;
end
clear t  

%data3= data2;
for j = 1: length(data)
for i = 11:2:(T-1)
    data2(j,i) = (data2(j,i-1)+data2(j,i+1))/2 ;
end
end

data= data2;
clear data2 
T=T-1;

%%%%%%%%
C = zeros(T,T);
Y = zeros(T,T);
D = zeros(T,T);

  % D Matrix: Indicator matrix if individual is present at both dates
  
  %diagonoal entries
  for j = 1:length(D)
  D(j,j) = sum( ~isnan(data(:,j+1)));
  end
  
  %off diagonal entries
  
  for j = 1:length(D)
      for i = j+1:length(D)
          resid = [data(:,j+1) data(:,i+1)];
            t=0;
            resid = resid';
            for l = 1:length(resid)
        if ~isnan(resid(:,l))
            t = t+1;
        end
            end
          D(j,i)=t;
         % D(j,i)=D(i,j);
          clear t resid
      end
  end
  
  % Y Matrix: Matrix of product of earnings of individual at both dates
  
     % diagonoal entries
    for j = 1:length(Y)
        help = data(:,j+1).*data(:,j+1);
        id1 = find(isnan(help));
        help([id1]) = [];
        Y(j,j) = sum(help);
        clear help id1
    end
    
     % off diagonoal entries

     for p = 2:T   
         for j = p+1:T+1
 resid = [data(:,p) data(:,j)];
 t=1;
    for i = 1:length(resid(:,1))
 if ~isnan(resid(i,:))
     resid2(t,:) =  resid(i,:) ;
    t = t+1;
 end
     end
 for l = 1: length(resid2(:,1))
 test(l) = resid2(l,2)*resid2(l,1) ;
 end
Y(p-1,j-1) = sum(test);
%Y(j-1,p-1) = Y(p-1,j-1);
clear resid resid2 test t 
end
end
  
% combine Y and D matrix

C= Y./D;


% stack C matrix

C=reshape(C,[(T)*(T),1]);  
      id1 = find(isnan(C));
        C([id1]) = [];
        
        
% MDE (with Identity matrix) 

% initial variance set to zero
var_z0 = (0)^2;

I = length(C);
   func = @(x) (C-theo(x(1),x(2),x(3),var_z0,0))'*eye(I,I)*(C-theo(x(1),x(2),x(3),var_z0,0)) ;
   x0 = [0.9,0.1, 0.1];
   lb = [0, 0, 0];
   ub = [.999, 5, 5];
   options = optimset('Display','iter','TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',10000);
   [x,fval,exitflag] = fminsearchbnd(func,x0,lb,ub,options);
   C_th = theo(x(1),x(2),x(3),var_z0,0);
   
   comp = [C,C_th];
   
   rho     = x(1);
   var_eta = x(3);
   var_e   = x(2);
   
   if gender == 11
       save('male_coup_labinc','rho', 'var_e', 'var_eta','var_z0','fval','C','C_th')
         elseif gender == 12
       save('female_coup_labinc','rho', 'var_e', 'var_eta','var_z0','fval','C','C_th')
        elseif gender == 13
       save('male_sing_labinc','rho', 'var_e', 'var_eta','var_z0','fval','C','C_th')
         elseif gender == 14
       save('female_sing_labinc','rho', 'var_e', 'var_eta','var_z0','fval','C','C_th')
   end
end
 



